Fractal clustering of inertial particles in random flows 



Jeremie Be<Q 

Institute for Advanced Study, Einstein Drive, Princeton, New Jersey 08540. 
Lab. G.-D. Cassini, Observatoire de la Cote d'Azur, BP 4^29, 06304 Nice Cedex 4, France. 

It is shown that preferential concentrations of inertial (finite-size) particle suspensions in turbulent 
flows follow from the dissipative nature of their dynamics. In phase space, particle trajectories 
converge toward a dynamical fractal attractor. Below a critical Stokes number (non-dimensional 
viscous friction time), the projection on position space is a dynamical fractal cluster; above this 
number, particles are space filling. Numerical simulations and semi-heuristic theory illustrating 
such effects are presented for a simple model of inertial particle dynamics. 



Dust particles, droplets, bubbles and various impurities 
advected by turbulent flow have usually a finite size and 
a mass density differing from that of the carrier fluid. 
Contrary to passive tracers, whose dynamics is conser- 
vative (when the flow is incompressible), the dynamics 
of such inertial particles is rendered dissipative by the 
Stokes drag, which can lead to strongly inhomogeneous 
spatial distributions. The full statistical description of 
such preferential concentrations is still an open question 
with many natural and industrial applications, such as 
the growth of rain drops in subtropical clouds the for- 
mation of planetesimals in the early Solar system^ op- 
timization of combustion processes and the coexistence 
problems between several species of plankton^ 

Maxey and Riley- derived an equation for the motion 
of a rigid spherical particle embedded in an incompress- 
ible flow. They assume that (i) the particle is smaller 
than the smallest turbulent scale of the carrier flow and 
that (ii) the Reynolds number associated to its size and 
its relative velocity with respect to the fluid is suffi- 
ciently small to approximate the surrounding flow by 
a Stokes flow. Then, the forces exerted on the parti- 
cle are buoyancy, the force due to the undisturbed flow, 
the Stokes viscous drag, the added mass effect and the 
Basset-Boussinesq history force. Because of the com- 
plexity of the resulting equation of motion, simpler mod- 
els are generally used. For instance, when the Stokes drag 
is very strong, the dynamics is close to that of passive 
tracer particles and the discrepancy can be captured by 
a spatial Taylor expansion, leading to a model in which 
the particles are advected by a synthetic flow compris- 
ing a small compressible component What singles out 
the model proposed here is its ability to take into ac- 
count the full phase-space dynamics of the particles and 
to capture the essential features of their dissipative mo- 
tion. We are interested in the "Batchelor regime" of the 
particles, meaning that we focus on spatial scales smaller 
than the Kolmogorov dissipation scale n. After rescaling 
of space, time and velocity respectively by factors n, rj 2 /v 
and v/rj, and assuming that the particle velocity is suf- 
ficiently close to that of the fluid, the Newton equation 
satisfied by its trajectory X(f) reduces to 
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where j3 = 3pf/(pf + 2p p ) is the added- mass factor (pf 
and p p are the fluid and the particle mass densities) and 
S, ; = a 2 /(3f3i] 2 ) is the Stokes number associated to the 
dissipative scale of the carrier flow (a being the particle 
radius). Introducing the co-velocity V = X — /3m(X, t), 
the equation of motion can be interpreted in terms of the 
(2 x (i)-dimensional dynamical system 



X = /?u(X,f)+V, 

V = ^[(l-/J)u(X,t)-V] 



(2) 
(3) 



The motion of the particles is clearly dissipative, even 
if the carrier flow is itself incompressible : indeed, when 
V • u = 0, the contraction rate in phase space reduces to 
—d/Sr) and is strictly negative, inducing a uniform con- 
traction. As a consequence, the long-time dynamics of 
the particles is characterized by the presence of an at- 
tractor, that is a dynamical fractal set of the phase space 
toward which the trajectories (X(t),V(£)) converge. Im- 
portant information on the dynamical system ©-Q, re- 
garding stability, Lyapunov exponents, etc is obtained 
from the linearized equation governing the separation 
R(i) = (<5X(t), 5*V(t)) between two infinitesimally close 
trajectories of the phase space. For scales within the vis- 
cous scale of turbulence, the velocity field u can be con- 
sidered spatially smooth and the separation R(t) obeys 
the linear differential equation 



R = 7W*R, M t = 



0<r(t) l d 
-5 — <r(t) —Id 
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where <x is the strain matrix of the carrier flow along the 
path of a reference particle: <Jij(t) = djUi(X(t),t), and 
Xd is the d-dimensional identity matrix. A full stability 
analysis of the dynamics can easily be doneji it relates 
the eigenvalues of the evolution matrix M. t to that of the 
stress tensor <x representative of the local structure of 
the carrier flow. In both two and three dimensions, this 
leads to distinguishing between heavy (/? < 1) and light 
particles (j3 > 1): the former are usually ejected from 
the elliptic regions, while the latter may cluster there in 
a pointwise manner. We therefore expect the vortex cores 
to be regions of high concentration of light particles and 
of low concentration of heavy particles, feature which is 
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generally observed in experiments and simulations (for a 
review, see Eaton and Fessler— ). 

We focus here on suspensions of particles with a volume 
fraction sufficiently small to neglect their interactions and 
collisions. Typically, the phase-space attractor on which 
the particles concentrate is a fractal object which may 
be characterized by various dimensions, in particular a 
non-random Hausdorff dimension du- As the position 
of the particles is obtained by projection from the 2d- 
dimensional phase space onto the d-dimensional physical 
space, the convergence to the attractor is responsible for 
strong inhomogcneities in the large-times distribution of 
particles. More precisely, a standard result of the geo- 
metrical theory of fractal sets^ states that if dn < d, 
the distribution of particles in the physical space is it- 
self a fractal set with Hausdorff dimension dn , whereas 
if dn > d, the particles fill the whole space. Hence, 
depending on the value of the dimension dn of the at- 
tractor, two different regimes are distinguished. Clearly, 
the dimension of the attractor is a function of the Stokes 
number and of the added mass parameter /3, and gen- 
erally also depends on the statistical properties of the 
velocity of the carrier fluid. Leaving aside this latter de- 
pendence, let us first note what can be easily inferred on 
the behavior of dn as a function of Stj and (3. First, in 
the limit of vanishing Stokes numbers, there is a reduc- 
tion of dimensionality and the dynamics of simple tracers 
is recovered. An initially uniform distribution of parti- 
cles remains uniform and we hence have dn — > d. Next, 
for very large Stokes numbers, the particles are less and 
less influenced by the carrier fluid and their motion be- 
comes ballistic. They thus fill the whole phase space and 
we have du — ► 2<i. In between these two asymptotic 
regimes, although it is not obvious that the dimension 
dn can fall below d, we shall actually show that there is 
a whole range of Stokes numbers for which du < d and, 
thus, preferential concentration on fractal clusters occurs 
in physical space. 

Finding theoretically or numerically the Hausdorff di- 
mension dn of the attractor is not, in general, a simple 
task: its determination demands a full understanding of 
the global dynamics and its numerical measurement re- 
quires very large numbers of particles. To obtain a simple 
estimate of the attractor dimension, Kaplan and YorkeiS 
proposed to use the Lyapunov dimension. It is given by 
the Lyapunov exponents Ai > . . . > \2d which measure 
the exponential growth of infinitesimal distances, sur- 
faces and volumes in the phase space and are expressed in 
terms of the limiting singular values of the Jacobi matrix, 
i.e. 

1 I* 
Xj = lim -\n\Jtej\ , J t = Texp / M s ds, (5) 
t^oo t J 

where the 's are the eigenvectors of the symmetric ma- 
trix Jt and Texp denotes the time-ordered exponen- 
tial. The Lyapunov dimension is defined as 
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where j satisfies Ai + . . . + Xj > and Ai + . . . + Aj+i < 
0. Beside being a simple estimate of the dimension of 
the attractor, it was actually shownii that the Lyapunov 
dimension gives a rigorous upper bound for the Hausdorff 
dimension dn of the attractor. 
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FIG. 1: Difference between the Lyapunov dimension d^ of heavy 
particles (/3 = 0) and the physical space dimension d, versus the 
Stokes number S,, (circle: d = 2, square: d = 3). The critical 
Stokes number corresponds to the value for which d^ = d. Upper- 
left inset : same in log-log coordinates showing a quadratic behavior 
at small Stokes numbers. 
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(a) S,, = 1CT 2 (b) S^ = 1 

FIG. 2: Snapshots of the position of N = 10 heavy particles 
(j9 = 0) associated to two different Stokes numbers: (a) smaller 
than the critical value, for which the particles form fractal clusters, 
and (b) larger than the critical value for which they fill the whole 
domain. The carrier incompressible flow is generated randomly by 
4 independent modes with a finite correlation time. 

We performed numerical experiments in two and three 
dimensions for a space-periodic carrier flow generated 
randomly by the superposition of few independent Gaus- 
sian Fourier modes with a correlation time of the order of 
unity (this specific form for the carrier flow was consid- 
ered by Sigurgeirsson and StuartjA2i who proved the exis- 
tence of a random dynamical attractor). The Lyapunov 
exponents are calculated by the use of the standard tech- 
nique of Benettin et ali&, and the resulting Lyapunov 
dimension is represented both for d — 2 and d = 3 in 
Fig. for the case of very heavy particles ((3 — 0). Two 
important observations can be made from these Simula- 
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tions. First, fractal clustering of particles already occurs 
at very small Stokes number where the Lyapunov dimen- 
sion behaves as g?l ~ d—CS^ with C > 0. This quadratic 
behavior near zero was predicted for the correlation di- 
mension by Balkovsky et al& using the method of the 
synthetic compressible flow cited earlier, as an approxi- 
mation at low Stokes numbers. The second observation 
is the presence of a critical value for the Stokes number 
below which the attractor dimension is smaller than d 
and where the particles form fractal clusters in the phys- 
ical space. The two regimes corresponding to different 
values of the Stokes number are illustrated for d = 2 in 
Fig. When the Stokes number is below the critical 
value (a), the particles concentrate onto a fractal set and 
both very dense and almost empty regions appear. On 
the contrary, for a Stokes number above the threshold 
(b) , the particles fill the whole domain, albeit with a non 
uniform density. 
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FIG. 3: Phase diagram in the parameter space S^) for the two- 
dimensional case, representing the three different regimes classified 
by the behavior of the Lyapunov exponents. 

When the particles have a finite mass (/3 ^ 0), there 
also exists a critical Stokes number for their concentra- 
tion onto fractal clusters. Figure shows for d = 2 the 
phase diagram obtained numerically which divides the 
parameter space {(3, S^) between three different regimes. 
When the sum of the two largest Lyapunov exponents 
is negative (light-gray domain), we have d^ < d and 
the particles form fractal clusters in the physical space. 
When the sum is positive (white domain), we have di > d 
and the particles fill the whole domain. The third case 
occurs only for particles lighter than the fluid and cor- 
responds to a negative largest Lyapunov exponent (dark 
gray area). The particles form pointwise clusters and, 
when the domain is bounded, they all converge to a sin- 
gle trajectory. 

We now present a heuristic argument, which we al- 
ready gave in the case j3 = of heavy particles?^ and 
which predicts the threshold in Stokes number. It relies 
on the use of the stability exponents, that are the ex- 
ponential growth rates of the eigenvalues of the Jacobi 
matrix J t defined in JSJ. Using Browne's theorem which 
bounds the singular values of a square matrix by its eigen- 



values, a necessary condition for the fractal clustering of 
particles is that the sum of the d largest stability expo- 
nents is positive. For heavy particles ((3 < 1), this sum 
can be estimated from the local analysis of the dynam- 
ics. First, since such particles tend to cluster within the 
hyperbolic regions of the flow, they are spending there a 
fraction of time much larger than in the elliptic regions. 
Let us assume that the relationship between the eigen- 
values of the evolution matrix Ait and those of the stress 
tensor cr(t) can be extended to the stability exponents, 
at least as an approximation; we can then derive a neces- 
sary condition for the presence of fractal clusters in the 
physical space. For d — 2, this condition can be easily 
written as 

S„ <^(/?-2 + 2v/l-/3 + /? 2 ), (7) 

where A/ is the (non-dimensional) largest Lyapunov ex- 
ponent associated to the carrier velocity field and calcu- 
lated along the trajectory of a simple fluid particle. It 
was easily verified that this bound is compatible with 
what is observed numerically in Fig. |3| 

It is often stated in the literature that the clustering 
of inertial particles is essentially due to the presence of 
coherent structures in turbulence (see., e.g., Squires and 
Eatoniii). It is indeed generally assumed that the struc- 
tures with long life times appearing in the flow are re- 
sponsible for a deterministic motion of the particle lead- 
ing to their concentration inside or outside the vortices. 
Although this argument, based on the local structure of 
the carrier flow, allowed us to find an upper bound for 
the critical value of the Stokes number, it is important 
to stress that preferential concentrations of particles arise 
solely as a consequence of the dissipative character of the 
motion. Indeed, we have also performed simulations with 
carrier flows which are delta-correlated in time and that 
are thus completely devoid of structured These simu- 
lations show that the dependence of the Lyapunov di- 
mension on the Stokes number S^ is very similar to that 
obtained in Fig. ^ The main difference is the behavior 
of g?£ as S, ; — » 0: in the delta-correlated case, the Lya- 
punov dimension tends linearly, and not quadratically, to 
the space dimension d. In the delta-correlated case, pref- 
erential concentrations are actually stronger than for a 
finite correlation time, contradicting the mechanism fre- 
quently invoked to explain concentrations. 

Of course, it is not sufficient to know that the parti- 
cles are concentrated in fractal objects; a fuller statistical 
description of their distribution is desirable. In particu- 
lar, for a quantitative description of their spatial inter- 
mittency, one needs their multifractal properties (scaling 
exponents of the various moments of the mass contained 
in a sphere of radius r; see e.g. Chap. 8 of the book by 
Frischi&). Preliminary results, for heavy particles, in- 
dicate that strong spatial intermittency can be present, 
even when the Lyapunov dimension is just slightly below 
that of the physical space. Let us also mention that more 
quantitative results are likely to be within reach using 
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multi-time methods in the asymptotics S n <C 1. Small 
Stokes numbers are of considerable interest since in most 
natural and industrial situations, this is the case. A last 
remark concerns the extension of this approach to the 
case of real turbulent flows. To confirm the existence of 
a threshold for fractal clustering of inertial particles, one 
needs to resolve scales which are much smaller than the 
Kolmogorov dissipation scale ry. This is quite a challenge, 
both for laboratory and numerical experiments. 

I am deeply grateful to K. Domelevo for initiating 
my work on inertial particles and to U. Frisch for his 



continuous support. During this work, I have benefited 
from stimulating discussions with A. Celani, M. Cencini, 
P. Horvai, K. Khanin, M. Stepanov and P. Tanga. Part 
of this research was possible thanks to the support of 
the Societe de Secours des Amis des Sciences which is 
warmly acknowledged. Part of this work was also sup- 
ported by the National Science Foundation under agree- 
ment No. DMS-9729992 and by the European Union un- 
der contract HPRN-CT-2000-00162. Numerical simula- 
tions were performed in the framework of the SIVAM 
project at the Observatoire de la Cote d'Azur. 



* Electronic address: bec@obs-nice. fr 

1 G. Falkovich, A. Fouxon and M.G. Stepanov, "Accelera- 
tion of rain initiation by cloud turbulence," Nature 419, 
151 (2002). 

2 A. Bracco, P.-H. Chavanis, A. Provenzale and E.A. Spiegel, 
"Particle aggregation in a turbulent Keplerian flow," Phys. 
Fluids 11, 2280 (1999). 

3 G. Karolyi, A. Pentek, I. Scheming, T. Tel and Z. 
Toroczkai, "Chaotic flow: the physics of species coexis- 
tence," Proc. Natl. Acad. Sci. USA 97, 13661 (2000). 

4 M.R. Maxey and J. Riley, "Equation of motion of a small 
rigid sphere in a nonuniform flow," Phys. Fluids 26, 883 
(1983). 

5 T. Elperin, N. Kleeorin and I. Rogachevskii, "Self- 
excitation of fluctuations of inertial particles concentration 
in turbulent fluid flow," Phys. Rev. Lett. 77, 5373 (1996). 

6 E. Balkovsky, G. Falkovich and A. Fouxon, "Intermittent 
distribution of inertial particles in turbulent flows," Phys. 
Rev. Lett. 86, 2790 (2001). 

7 J. Bee, Particules, singularites et turbulence (PhD thesis, 
Universite P. & M. Curie, Paris, 2002), also available at 
http: / /www.obs-nice.fr /bee /these. pdf 

8 J.K. Eaton and JR. Fessler, "Preferential concentrations 
of particles by turbulence," Int. J. Multiphase Flow 20, 
169 (1994). 



9 K.J. Falconer, The geometry of fractal sets (Cambridge 
University Press, Cambridge, 1986). 

J.L. Kaplan and J. A. Yorke, "Chaotic behavior of multidi- 
mensional difference equations," in Proc. Functional differ- 
ential equations and approximation of fixed points (Bonn, 
1978), Lecture Notes in Math. 730 (Springer, Berlin, 1979) 
pp. 204-227. 

1 A. Douady and J. Oesterle, "Dimension de Hausdorff des 
attracteurs," C. R. Acad. Sc. Paris 24, 1135 (1980). 

2 H. Sigurgeirsson and A.M. Stuart, "A model for preferen- 
tial concentration," Phys. Fluids 14, 4352 (2002). 

3 G. Benettin, L. Galgani, A. Giorgilli and J.M. Strelcyn, 
"Lyapunov characteristic exponents for smooth dynamical 
systems and for Hamiltonian systems," Meccanica 15, 9 
(1980). 

4 J. Bee and K. Domelevo, "Inertial particles in random 
flows," in Advances in Turbulence IX, Proc. Ninth Eu- 
ropean Turbulence Conference (Southampton 2002), LP. 
Castro and RE. Hancock Eds. (CIMNE, Barcelona, 2002) 
pp. 31-34. 

5 KD. Squires and J.K. Eaton, "Preferential concentration 
of particles by turbulence," Phys. Fluids A 3, 1169 (1991). 

6 U. Frisch, Turbulence: the Legacy of A.N. Kolmogorov 
(Cambridge University Press, Cambridge, 1995). 



